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Abstract 

The Allen-Cahn equation is solved numerically by operator splitting Fourier spec¬ 
tral methods. The basic idea of the operator splitting method is to decompose the 
original problem into sub-equations and compose the approximate solution of the 
original equation using the solutions of the subproblems. Unlike the first and the 
second order methods, each of the heat and the free-energy evolution operators has 
at least one backward evaluation in higher order methods. We investigate the effect 
of negative time steps on a general form of third order schemes and suggest three 
third order methods for better stability and accuracy. Two fourth order methods are 
also presented. The traveling wave solution and a spinodal decomposition problem 
are used to demonstrate numerical properties and the order of convergence of the 
proposed methods. 

Key words: Operator splitting method; Allen-Cahn equation; Heat evolution 
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1 Introduction 


The Allen-Cahn (AC) equation was originally introduced as a phenomenolog¬ 
ical model for antiphase domain coarsening in a binary alloy [1]: 


dt 


F'{<P{yiV)) 

^2 


+ A0(x,t), 


X G 0 < t < T, 
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where O is a domain in {d = 1,2,3). The quantity ^(x,t) is defined as 
the difference between the concentrations of two components in a mixture, for 
example, 0(x, t) = (m^ — rnis)/(ma + m-is) where and are the masses 
of phases a and /5. The function F(0) = 0.25(0^ — 1)^ is the Helmholtz free- 
energy density for which has a double-well form, and e > 0 is the gradient 
energy coefficient. The system is completed by taking an initial condition 
0(x, 0) = 0°(x) and a homogeneous Neumann boundary condition V0-n = 0, 
where n is normal to dVt. 

The AC equation and its various modified forms have been applied in ad¬ 
dressing a range of problems, such as phase transitions [1], image analysis 
[2,3], motion by mean curvature [4,5,6], two-phase fluid flows [7], and crystal 
growth [8,9,10]. Therefore, many researchers have studied numerical meth¬ 
ods for solving the AC type equation to improve stability and accuracy and 
to have a better understanding of its dynamics. Stable time step size of ex¬ 
plicit schemes is severely restricted due to the nonlinear term F'{(p) and im¬ 
plicit schemes suffer from a solvability problem with large time steps. One of 
considerable semi-implicit methods is unconditionally gradient stable method 
proposed by Eyre [11], which is first order accurate in time, and uncondition¬ 
ally gradient stable means that a discrete energy non-increases from one time 
level to the next regardless of the time step size. And the authors in [12,13] 
proposed first and second order stabilized semi-implicit methods. 

Another numerical method employed for solving the AC equation is the oper¬ 
ator splitting method [12,14,15]. Operator splitting schemes have been applied 
for many types of evolution equations [16,17,18,19,20,21]. The basic idea of 
the operator splitting method is to decompose the original problem into sub¬ 
problems which are simpler than the original problem and then to compose 
the approximate solution of the original problem by using the exact or ap¬ 
proximate solutions of the subproblems in a given sequential order. Operator 
splitting methods are simple to implement and computationally efficient to 
achieve higher order accuracy while semi-implicit schemes are hard to im¬ 
prove the order of convergence. The first and the second order operator split¬ 
ting methods for the AC equation is quite well-known [12,14,15], however, the 
higher order (more than two) operator splitting method for the AC equation 
is less well-known. 

In this paper, we investigate higher order operator splitting schemes and pro¬ 
pose several higher order methods to solve the AC equation with a Fourier 
spectral method. We decompose the AC equation into heat and free-energy 
evolution equations, which have closed-form solutions in the Fourier and phys¬ 
ical spaces, respectively. Because the first and second operator splitting meth¬ 
ods have only forward time steps, the boundedness of the solution is guar¬ 
anteed regardless of the time step size [15]. However, we could not guarantee 
the stability with large time step size since each operator has at least one 
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backward time step with third and higher order of accuracy [17,18]. Because a 
backward time marching affects numerical stability on both sub-equations, we 
consider ways of minimizing the effect of negative time steps and introduce a 
cut-off function to limit the exponential amplification of high-frequency modes 
in solving the heat evolution equation. 

This paper is organized as follows. In section 2, we briefly review the operator 
splitting methods which are studied by the authors in [17]. In section 3, we 
present higher order operator splitting Fourier spectral methods for solving the 
AC equation. We discuss the stability issues for backward time marching and 
suggest the three third order operator splitting methods. We present numerical 
experiments demonstrating numerical properties and the order of convergence 
of the proposed methods in section 4. Conclusions are drawn in section 5. 


2 A brief review on the operator splitting method 


In this section, we review some of the basic properties of the operator splitting 
methods for a time evolution equation with two evolution terms in summa¬ 
rizing the work by D. Goldman and T. Kapper [17]. Let be the solution 
operator for the time evolution equation = /a(^), that is := 

-|- a At), and be the solution operator for fB{<f>)- Then the operators 
A and B satisfy the semi-group properties. Suppose we want to minimize the 
number of the operator evaluations of and B^^^ in order to get a N-th 
order approximation of the following ordinary differential equation consists of 
two evolution terms, 

^ = nw+/ b(«. (2) 

It is well-known that the simplest form of the first order solution operator for 
(2) is given as 

5 ( 1 ) ^ ^ 3 ) 

that is, is ^ first order accurate approximation of 0(t-|-At). Here the 

choice of A and B (or fA and fs) is arbitrary, thus without loss of generality, 
we may assume that the first operator evaluated is always 


We now consider a solution operator with 2p (or 2p—1 if bp = 0) evalua¬ 
tions of the operators A and B in the following form, 

^(p) — J^bpAt j^apAt _ _ _ ^biAt ^aiAt 

where all of {hj}^Z^ are non-zeros. The coefficients ai,...,ap and 

bi,... ,bp in must satisfy certain conditions to make an A-th order 
approximation operator for (2). It is well-known that there exists at least 
A-th order accurate when p > N. (See [17] and the references therein for 


3 



the derivation of the following conditions.) For first-order accuracy, {aj}, {bj} 
must satisfy 

p p 

= = ( 5 ) 

j=i j=i 

For second-order accuracy, {aj} and {bj} must satisfy (5) and the conditions 


P \ P ( ^ \ 1 

=^- (6) 

i =2 \k=l / j=l \k=i / ^ 

For third-order accuracy, {aj} and {bj} must satisfy (5), ( 6 ), and the condi¬ 
tions 

(7) 


^J-1 


— X] bj 


j =2 \k=l 


J=1 


For a second-order scheme of the form (4) with p = 2, ^a 2 At At 

^aiAt, g-.^g 

Ui -|- n2 = 1, bi -{- b2 = 1, n2^i = 2’ (®) 

Since there are three equations for the four unknowns, let bi=u) 0 ) be a 
free parameter, then the solution of ( 8 ) gives a general form of a second order 
solution operator with up to 4 operator evaluations. 

Note that A^ with cu = 1 is the simplest form (with only 

three evaluations) among second order operators since two evaluations of A 
and B is not enough to make it second order accurate. 


For a third-order scheme of the form 


5(3) 


_ jgbsAt ^asAt 2^62 At ^a 2 At j^hiAt ^aiAt 


( 10 ) 


(5), ( 6 ), and (7) give 


ai + a 2 + as — 1, bi + b 2 + bs — 1, 0261 + 03(^1 + ^2) — 

02^1 + 03(^1 + ^2)^ = g) bia{ + 62(01 -|- 02)^ + 63 = -. (12) 

Choosing 63 = cu to be a free parameter, we can obtain two branches of the 
solution for ( 11 ) and ( 12 ), 


1 —uj \J D{u) 4n; — 1 ^ | — 6f 02 

^ 2 ^ 2(4n; - 1) ’ ^ 2(3n; - 1) ’ ^ 1 - uj 
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where 


bf = 1 - bf - 63 , 


af = 1 — 02 — Og , 


D{u) = {cu — 1 Y{Auj — lY + 12(4a; — 1) (a; 



Note that real solutions of (11) and (12) are only possible for cn > ^ and 
uj < u)*, where oj* ~ —1.217 • • • is the real root of D{u))/(Aoi) — 1) = 0. 
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Fig. 1. Positive branch solutions, 6 ]'', a^, 6 ^, 6 ^ as a function of 6 ^ — lu 

Figure 1 shows the positive branch solutions, a]*", 61 ", a^, 6 ^, 03,63 as a func¬ 
tion of 63 = LJ for the third order operator <S^+ and Figure 2 shows the 
negative branch solutions for S^-. In any case, there exists exactly one neg¬ 
ative value among ai, a 2 , ag and also only one negative value among 61 , 62 , & 3 - 
There are three special cases when the solutions may blow up. As a; — 

with 02 = 0 , 62 + = I degenerates into a second order operator, 

A^^*. As (u —> |, S^l with b^ = 0, af + \ or b^ = 0, 

(I 2 A-ds =1 degenerates into a second order operator, A^^^ Bs^^ A^^^. 
As a; 1 , the negative branch solutions have removable singularities and S^- 
converges to B^^ A'^^* B^^* A^^^ Bi^^ A^^^^ whereas the positive branch 
solution does not provide a convergent operator. 

We remark that a symmetric with 63 — 0, ai = ag, and 61 = 62 satisfying 
(5), ( 6 ) has only a second-order accuracy, that is, Ui = |, 61 = |, and ^2 = 1 
does not satisfy (7). However, a symmetric with 64 — 0, Otl — 0 . 4 , 0,2 — 0-3; 
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Fig. 2. Negative branch solutions, , 6 ^^ , a 2 , 62 ? <^3 5 ^3 a function of 63 — uj 
and hi = 63 satisfying only (5), ( 6 ), and (7), 

5 (f) := g{l-20.)At j^coAt At (-^ 3 ^ 

happens to be a fourth-order accuracy with cu — uu=l/{‘2 — 2^/^) 1.3512, 

—0.1756, and 1 — 2a; ~ —1.7024. This is the simplest form of fourth 
order operator with only 7 operator evaluations and this can be derived as a 
symmetric combination of a second order operator 7”^* A~ A~, 

:= 7-(l-2a;)At ^uAt^ UJ = UJ^. 

Another a well-known fourth order operator splitting method [22] can be also 
derived as a symmetric combination of the second order operator 7”^*, 

^(6)_/y-'CdAt /y-'CdAt 4c(;)At 'j^uoAt ^CdAt 

_ At j^uo At At j^uo At At'j^{l—Auj)At At At At j^uo At At 

with UJ = a;y=l/(4 — 4^/^) ~ 0.4145. The (Sy ^ method is computationally less 
efficient (11 operator evaluations compared to minimum of 7 evaluations) but 
has better stability condition ~ —0.1217, 1 — 4a; ~ —0.6580) than the 
method defined in (13). 

We close this section with a remark that not just the third and the fourth 
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order methods mentioned above but any operator splitting methods of third or 
higher order contains at least one negative time steps for each of the operators, 
A, B. (See [17,18] for the proof of the general theorem.) 


3 Higher-order operator splitting Fourier spectral methods 


We consider the AC equation (1) in one-dimensional space fl = {0,L). Two- 
and three-dimensional spaces can be analogously defined. For simplicity of 
notation, we sometimes abuse the notation 0 = (pit) referring 0(-, t) and define 
the free-energy evolution operator’^ as follows 

W^*(0(r)):=0(r + At), (15) 

where 0(t” -|- At) is a solution of the first order differential equation 

^ _ F'{(P) 

dt A 

with an initial condition 0(t”). For given F'icp) = <jA — <p, we have an analytical 
formula (See [12,14,15]) for the evolution operator in the physical space 




_ f _ 

\/(/)^ -h (1 - (lF)e~~^ 


(16) 


We also define the “heat evolution operatoF as follows 


n^\(PiF)) := (/)(r + At), 


(17) 


where + At) is a solution of the first order differential equation 



with an initial condition In this paper, we employ the discrete cosine 

transform [23] to solve the AC equation with the zero Neumann boundary 
condition: for A; = 0,, M—1, 


M-l 

= OLk 4^1 COS 
/=0 


^.k(l+^ 
M 




where (pi = (p and ag = for k 

have a semi-analytical formula for the evolution operator 77^* 
cosine space 


n^\(P) = c-^ 


[f] 


> 1. Then, we 
in the discrete 

(18) 
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where Ak 



and C denotes the discrete cosine transform. 


For the first order operator splitting scheme in (3) and the second order 
scheme in (9) with 0 < a; < 1, the evaluations are all forward time 
marching, that is, all of and are positive. We can easily show 

that both schemes are unconditionally stable, in the sense that |0(F^+At)| < 1 
if |0(t”)| < 1 regardless of the time step size. (See [15] for the proof.) However, 
in the case of third or higher order, each of operators T, T-L has at least one 
backward evaluation as mentioned in section 2. For this reason, we need to 
investigate the stability of the operators and especially for large 

At. 


The stability issue for backward time heat equation is well-known. Even though 
-^ibAt -^=FAt (^^ithout noise) is always an identity operator regardless of the 
size of At, the numerical composition of the operators (even with small error) 
is far away from the identity operator when At becomes large since is 

exponentially big for At ^ 1. The stability of the numerical composition of 
the free energy evolution operators is less well-known and we want to explain 
why the numerical composition of the operators (even with small 

error) is far away from the identity operator when At becomes large using the 
following figure. 



Fig. 3. </>(At) = ^^*(^(0)) with various initial values, ^(0) = —1.6, —1.4, • • • , 1.6 

Figure 3 plots J^^*(0) as a function of At with various initial values of (p 
between —1.6 to 1.6. As you can see, W^*(0) with |0| < 1 converges to ±1 as 
At 1, however, the solution with |0| > 1 as a result of small perturbation 
may blow up when At <C —1. Therefore, composition of two operators 
followed by even with small evaluation error near 1 is no longer bounded 
as At is getting bigger. And J^“^*((/)) with |(/)| < 1 converges to 0 for At 1 
thus followed by for At 1 may change the sign of result even 

with small perturbation near 0. This non-linear stability effect is basically a 
consequence of the fact that the solution of the free-energy evolution operator 
is exponentially close to 1 or 0 as At ^ ±oo. 
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Fig. 4. Minimum and maximum of {a^and {b '^as a function of 6 ^ = u. 
The region where values are bounded by [—1,1] is shaded in yellow. 



Fig. 5. Minimum and maximum of {a^ and {bj as a function of 63 = lo. 
The region where values are bounded by [—1,1] is shaded in yellow. 


In order to achieve a better stability condition, we propose third order schemes 
with bounded values of {aj,bj}j^i. Figures 4 and 5 show the minimum and 
maximum values of and for the positive and negative branches, 

respectively. It is worth noting that af ,bi ,a 2 ,b 2 ,a^, 63 are bounded by [—1,1] 
when 0.26376-•• < < 0.29167 •• • and af, 67 , > ^3 ^-^e bounded 

by [—1,1] when 0.26376 ■ ■ ■ < uj~ < 0.27362 • • • or 1/2 < < 1. Since there 

is exactly one negative value among max{oj}|^^ > — min{aj}|^^ can 

be inferred from (5) when \aj\ < 1. Similarly max{ 6 _,}|^;^ > —mm{bj}^^i 
when \bj\ < 1. In the shaded regions on the hgures where values of \aj\, \bj\ 
are bounded by 1, there are three local minima of max{|aj|, at which 

values are summarized on the following table. 


Table 1 

Solutions for , a^,b^,a^,b^,af,bf 


at the local minima of max{ 





Condition 

ai 

hi 

a2 

b2 

as 

hs 

(jJX 

a+ = fe+ 

0.78868.. 

-0.07189.. 

-0.44191.. 

0.78868.. 

0.65324.. 

0.28322.. 

UJY 

fef = aj 

0.26833.. 

0.91966.. 

-0.18799.. 

-0.18799.. 

0.91966.. 

0.26833.. 

UJZ 

®2 “ ^3 

0.28322.. 

0.65324.. 

0.78868.. 

-0.44191.. 

-0.07189.. 

0.78868.. 


It is worth to note that the sets of {aj and {bj are same for oj = ujy. 
The set {adfor uj^ = u)x is {bj for u)~ = u)z and the set (6+ for 
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u)^ = ojx is {aj for cu~ = uz- This symmetry gives us a freedom to choose 
the order of operator evaluations and we define three third order operator 
splitting methods Sx\Sy\Sz^ for the AC equation as follows: 


^(3) _ '^asAt j^b2At '^a2At j^biAt ^aiAt 


(19) 


where and {bj}^^i are given in Table 1. 

Another issue raised with negative time step is that the heat evolution op¬ 
erator < 0 may amplify the high frequency modes exponentially 

big, S> 1. This situation —A^At = (jA'j At ^ 1 happens when 

k^At S> 0(1). On the other hand, a physically reasonable bound for At in the 
AC equation is ^ < 0(1), thus the blow-up may occur only for physically too 
high frequency modes, A: S> Thus, we introduce a cut-off function to bound 
of for high frequency modes where —A^At ^ 1 . We will numerically 

demonstrate the effect of introducing the cut-off function in subsection 4.1. 


4 Numerical experiments 


In this section, we numerically demonstrate the order of convergence of the 
proposed third order schemes in (19) and the fourth order 

schemes <5^^ in (13) and in (14). Two examples are used for the test, 
one is the traveling wave solution with analytic solution and the other is a 
three-dimensional spinodal decomposition problem with random initial val¬ 
ues. 


One of the well-known traveling wave solutions of the Allen-Cahn equation is 


^{xA) 


1 

2 


^1 — tanh 


X — 0.5 — st\ 

2V2e ) ’ 


( 20 ) 


where s = 3/(\/2e) is the speed of the traveling wave. The leftmost plot in 
Figure 6 shows the initial profile 0(x, 0) and the analytic solution (/){x, Tf) at 
Tf = 1/s with e = 0.03^/2. Using this traveling wave solution, we compare 
the first, second, third, and fourth order operator splitting Fourier spectral 
methods described in section 3. The numerical solutions (f){x,t), 0 < t < Tf 
are obtained with various time step sizes At but the spatial grid size is fixed to 
h = 2~^ which provides enough spatial accuracy. The traveling wave solution 
with the same numerical parameters are used in the following two subsections 
to test the third and the fourth order schemes. 


The rightmost plot in Figure 6 shows the numerical error of the first order 
scheme in (3) and the second order scheme in (9) compared to the 
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Fig. 6. Traveling wave solution (f){x,Tf) at Tf — Ijs with e = 0.03\/2. And relative 
I 2 errors of 4>{x,Tf) by with h = 2~^ for various time step sizes At. 

analytic solution at t = T/. It is worth to remind that the first and the second 
order schemes apply only forward time steps of T and H, thus the stability 
(or boundedness of the solution) regardless of the size At can be easily proven. 
(See our previous paper [15] for numerical properties of these first and second 
order schemes.) 


4.1 Cut-off function and stability of the third order methods 


As mentioned in section 3, negative time steps of 4F and l-L are unavoidable 
in the third or higher order operator splitting methods. Especially a negative 
time step makes the heat evolution operator exponentially big, therefore, we 
introduce the following cut-off function with a tolerance Ktoi for the heat 
evolution operator "H, 


= c-^ 


Cif] 


( 21 ) 


The choice of Ktoi depends on the time step size OjAt and highest frequency 
modes k^ax which are functions of desired computational accuracy. Following 
computational examples in this subsection give a basic guideline for the choice 
of Ktoi- 


As mentioned in section 2, we have various coefficients and as 

a function of = m. To investigate the effect of u) in the third order method 
51 ?' , we consider the traveling wave problem given in (20). We compute relative 
I 2 errors for various u values with a fixed time step At = 2~^/s and Figures 7 
(a) and (b) show relative I2 errors of the traveling wave solution 0(x, Tf) by 
the third order methods for positive and negative branches of various cn, 
respectively. Here we set Ktoi = 10^ (blue solid fine) or 10® (green dashed fine). 


The first noticeable point in Figure 7 might be that the error is relatively 
large at or ^ | where the third order operator degenerates into 
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(a) With positive branch solutions for in (10) 



Ktoi = 10 " 
o- - Ktoi =10^ 

Snnnoeeeeo 





-2 - 1.8 - 1.6 - 1.4 - 1.2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

-2 <LJ <LJ* 1/4 < a; < 2 

(b) With negative branch solutions for in (10) 


Fig. 7. Relative I 2 errors of the traveling wave solution (f){x,Tf = 1/s) by the third 
order method for various u with At = 2“^/s, e = 0.03\/2, and h = 2~^. 


a second order operator. Also a region near = 1 in the positive branch 
case, the computation does not provide any accuracy at all. As 1, S'^'^ 

contains a big negative time step of the heat evolution operator since 

min{aj} —>■ — 00 . In these cases, the choice of cut-off parameter Ktoi becomes 

important, and small Ktoi is recommended when — min{oj}At ::§> ^ . 

For 0.26376 • • • < < 0.29167 • • • in which {a/} and {6/} are bounded 

by [—1,1], especially near ux at which max{|oj|, \bj\} has a local minimum, 
the error is smaller than that for other u values. The similar phenomenon is 
observed the computation for the negative branch. We choose three special 
values = cox, <->J~ = uJy, and (jJ~ = ujz for Sx\ <5y/ and Sz \ respectively. 
For these cases, all {aj} are bounded by [—1,1] and the choice of cut-off value 
Ktoi does not play an important role in the computation. 


We now investigate the effect of highest frequency k^ax to Ktoi- Plots in Fig¬ 
ure 8 show relative I 2 errors of the traveling wave solution 0(x, Tf) by the 
third order method ^ with different spatial grid sizes h = ^ = ^ = 2“® 
or h = = 2“®. If a cut-off function is not used (labeled as Ktoi = Inf), 

the computation provides no accuracy for relatively large time step. The com¬ 
putation may even stop as two biggest At cases for M = 1024 and the cases 
happen more often as kmax = M becomes large. If At is larger than e^, Ktoi 
must be properly chosen in order to valence the accuracy loss while avoiding 
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blow-up. However, the choice of K^oi makes no significant difference of the 
solution when At < (which is physically valid limit for the AC equation) 
since the high frequency modes with ^ is negligible for the physically 
meaningful solution of the AC equation. So the simplest rule of thumb might 
be setting Ktoi around the desired accuracy of the computation. 


(a) M = 256 



Fig. 8. Relative I 2 errors of the traveling wave solution 4>{x,Tf = 1/s) by Sy with 
e = 0.03^2. 



Time step 


-10 
^ 10 ' 
= 10 ^' 
- Inf 


4-2 Convergence of the third and the fourth order methods 


We implement the proposed third order schemes (Sy \ (S^ ^ in (19) and the 
fourth order schemes Sjj^ in (13) and (Sy^ in (14). We set the spectral grid 
size h = 2“®, the cut-off limit Ktoi = 10® and compare the numerical solutions 
for various time step sizes At with the analytic solution (20) with e = 0.03v^. 
Figure 9 numerically indicates that the proposed methods have the third and 
the fourth order accuracy, respectively. 




(a) (b) 


Fig. 9. Relative I 2 errors of the traveling wave solution (j){x^Tf) at Tj = I /5 by (a) 
the third order methods (b) the fourth order schemes S^\ 
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4-3 Convergence of the spinodal decomposition problem in 3D 


In this subsection, we compute a spinodal decomposition problem satisfying 
the AC equation (1) in three-dimensional space with e = 0.015. The intervals 
(—1,—1/\/3) and (l/-\/3,1) where > 0 are called metastable intervals 

and (—1/\/3, l/\/3) where F"((f)) < 0 is called the spinodal interval [24]. It is 
known that f which lies in the spinodal interval is very unstable and the growth 
of instabilities results in phase separation, which is called spinodal decompo¬ 
sition. In order to check the numerical convergence, we integrate (f){x, y, z, t) 
up to time Tf = 0.01 by the proposed numerical schemes with various time 
step sizes At = 10“^/2, • • • , 10“^/2^. The initial condition is given on the 
computational grid with h = 2“® in the domain Q = [0,1] x [0,1] X [0,1] 
as = 0.005 • rand(x,y,z) where rand(x,y,z) is a random number 

between —1 and 1. Figure 10 shows the initial and the reference solutions at 
t = 10“^, 10“^ computed by the fourth order numerical scheme ^ with the 
numerical parameters Ktoi = 10^ and At = 10“^/2^. 



t = 0 t = 10-^ t = 10-2 

Fig. 10. The reference solutions (l){x^y^ z^t) by the fourth order method with 
Ktoi = 10^ and At = 10-^2^ 






10“® 10"® 10"“ lO""^ 

Time step size 

(b) Ktoi = 10'’ 


Fig. 11. Relative I 2 errors of z^Tf = 0.01) by ^ 

^ with various time step sizes At = 10-^/2, • • • , 10-^/2^. 
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We also implement the first order scheme in (3), the second order scheme 
in (9), the proposed third order schemes \ <Sy \ ^ in (19), and the 

fourth order schemes in (13) and in (14). The numerical results in 
Figure 11 show that the cut-off value Ktoi does not play a role when At 
is smaller than while the computational results have marginal difference 
when At is greater than e^. The accuracy results numerically demonstrate the 
proposed schemes provide the expected order of convergence in time. 


5 Conclusions 


We proposed and studied the higher order operator splitting Fourier spec¬ 
tral methods for solving the AC equation. The methods decompose the AC 
equation into the subequations with the heat and the free-energy evolution 
terms. Unlike the first and the second order methods, each of the heat and 
the free-energy evolution operators has at least one backward evaluation in 
the higher order methods. For the third order method, we suggested the three 
values ux,0i>Y,‘^z af which max{|aj|, |6j|} have local minimums and we then 
obtained smaller error than other uj values. For the fourth order method, we 
used two symmetric combinations of the second order operators. And a sim¬ 
ple cut-off function could limit exponential amplification of the high frequency 
modes in the heat operator and it worked well with the proposed schemes. We 
numerically demonstrated, using the traveling wave solution and the spinodal 
decomposition problem with random initial values, that the proposed methods 
have the third and the fourth order convergence as expected. 
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